Importing GPS data into R
require(raster);require(sp);require(rgdal); require(ggmap)
## Loading required package: raster
## Loading required package: sp
## Loading required package: rgdal
## rgdal: version: 1.2-18, (SVN revision 718)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.3, released 2015/09/16
## Path to GDAL shared files: /usr/share/gdal/1.11
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-7
## Loading required package: ggmap
## Loading required package: ggplot2
# Set your working directory
setwd('/home/pjg/GIS')
# Create a directory folder and move into it
dir.create('AdvancedGIS_R')
## Warning in dir.create("AdvancedGIS_R"): 'AdvancedGIS_R' already exists
setwd('AdvancedGIS_R')
If you are using data from a Garmin GPS
##############################################################
############### Importing Garmin GPS data ##################
##############################################################
###### Create a character vector of the entire pathway for the gpx object from your Garmin GPS
gpfile<-"/home/pjg/GIS/Classes/spatial_bioinformatics/R_GIS/Current.gpx"
###### Read in the GPX file using the function from rdgal
GPSfile<-readOGR(dsn = gpfile, layer = "track_points")
## OGR data source with driver: GPX
## Source: "/home/pjg/GIS/Classes/spatial_bioinformatics/R_GIS/Current.gpx", layer: "track_points"
## with 3 features
## It has 26 fields
##### Extract all of the data and coordinates
wpsDF<-GPSfile@data
wpsCoords<-GPSfile@coords
##### Pull out only the necessary data:longitude, latitude, time, elevation
##### Garmin does not supply accuracy at each waypoint
wps<-cbind(wpsCoords[,1:2], wpsDF[,c(5,4)])
colnames(wps)<-c("lon", "lat", "time", "elev")
head(wps)
## lon lat time elev
## 1 -73.94220 40.85085 2018/05/06 20:09:54+00 81.85
## 2 -73.94218 40.85087 2018/05/06 20:10:28+00 81.85
## 3 -73.94220 40.85085 2018/05/06 20:10:57+00 81.37
If you are using GPS data from Cellphone
###############################################################
############ Importing cellphone data #######################
###############################################################
# Read in .csv file of waypoints
wps<-read.csv('/home/pjg/GIS/Classes/spatial_bioinformatics/R_GIS/20180413.csv')
# Take a look at the first few rows of the data.frame
head(wps)
## time lat lon elevation accuracy bearing
## 1 2018-04-13T15:01:31.966Z 41.22050 -74.29108 NA 20.274 NA
## 2 2018-04-13T15:03:20.109Z 41.22051 -74.29101 NA 20.249 NA
## 3 2018-04-13T15:03:20.110Z 41.22051 -74.29101 NA 20.249 NA
## 4 2018-04-13T15:04:47.000Z 41.22035 -74.29116 263 37.000 NA
## 5 2018-04-13T15:05:08.000Z 41.22019 -74.29070 254 39.000 NA
## 6 2018-04-13T15:05:36.119Z 41.22050 -74.29095 NA 20.460 NA
## speed satellites provider hdop vdop pdop geoidheight ageofdgpsdata
## 1 NA 0 network NA NA NA NA NA
## 2 NA 0 network NA NA NA NA NA
## 3 NA 0 network NA NA NA NA NA
## 4 NA 11 gps 1.2 0.9 1.5 -37 NA
## 5 0 11 gps 1.2 0.9 1.5 -37 NA
## 6 NA 0 network NA NA NA NA NA
## dgpsid activity battery annotation
## 1 NA NA 90 NA
## 2 NA NA 90 NA
## 3 NA NA 90 NA
## 4 NA NA 90 NA
## 5 NA NA 90 NA
## 6 NA NA 89 NA
# We are only really interested in the first 5 columns of the waypoint file, but let's rearrainge them a bit
wps<-wps[c(3,2,1,4,5)]
# Now lon and lat are the first 2 columns of the data.frame
head(wps)
## lon lat time elevation accuracy
## 1 -74.29108 41.22050 2018-04-13T15:01:31.966Z NA 20.274
## 2 -74.29101 41.22051 2018-04-13T15:03:20.109Z NA 20.249
## 3 -74.29101 41.22051 2018-04-13T15:03:20.110Z NA 20.249
## 4 -74.29116 41.22035 2018-04-13T15:04:47.000Z 263 37.000
## 5 -74.29070 41.22019 2018-04-13T15:05:08.000Z 254 39.000
## 6 -74.29095 41.22050 2018-04-13T15:05:36.119Z NA 20.460
Visualizing waypoints on satellite imagery
###############################################################
############ Visualizing and creating features ##############
###############################################################
## Visualizing your waypoints with satellite data
# Find the bounding box of your waypoints
loc = cbind(c(min(wps$lon), max(wps$lon)), c(min(wps$lat), max(wps$lat)))
# Name the columns
loc = as.data.frame(loc)
colnames(loc) = c('lon', 'lat')
# Get the images from Google Earth Engine
manbox <- make_bbox(lon = loc$lon, lat = loc$lat, f = .1)
manmap <- get_map(location = manbox, maptype = "satellite", source = "google", zoom =18)
# Plot
ggmap(manmap) +
geom_point(data = loc,
color = "red",
size =1) # for better resolution, change size to smaller
Creating Shapefile Features from localities
## Transform the waypoints and save as a shapefile
wpsShp<-SpatialPointsDataFrame(coords = wps[,1:2], data = wps)
#writeOGR(obj = wpsShp, dsn = paste(getwd()), layer = "Allwaypoints", driver = "ESRI Shapefile")
## Separate the data into the appropriate features from which waypoints were taken outside
The first step is to take another look at the data and find the rownumber of the corresponding feature
print(wps)
## lon lat time elevation accuracy
## 1 -74.29108 41.22050 2018-04-13T15:01:31.966Z NA 20.274
## 2 -74.29101 41.22051 2018-04-13T15:03:20.109Z NA 20.249
## 3 -74.29101 41.22051 2018-04-13T15:03:20.110Z NA 20.249
## 4 -74.29116 41.22035 2018-04-13T15:04:47.000Z 263 37.000
## 5 -74.29070 41.22019 2018-04-13T15:05:08.000Z 254 39.000
## 6 -74.29095 41.22050 2018-04-13T15:05:36.119Z NA 20.460
## 7 -74.29109 41.22045 2018-04-13T15:06:36.470Z NA 23.969
## 8 -74.29080 41.22039 2018-04-13T15:07:19.000Z 216 38.000
## 9 -74.29077 41.22048 2018-04-13T15:06:54.837Z NA 21.048
## 10 -74.29083 41.22043 2018-04-13T15:07:36.000Z 204 37.000
## 11 -74.29064 41.22059 2018-04-13T19:06:07.000Z 180 21.000
## 12 -74.29064 41.22055 2018-04-13T19:07:38.000Z 186 40.000
## 13 -74.29066 41.22055 2018-04-13T19:08:47.000Z 188 34.000
## 14 -74.29067 41.22056 2018-04-13T19:10:09.000Z 186 31.000
## 15 -74.29068 41.22057 2018-04-13T19:10:13.000Z 185 16.000
## 16 -74.29068 41.22057 2018-04-13T19:10:16.000Z 185 18.000
## 17 -74.29120 41.22047 2018-04-13T19:09:43.811Z NA 19.906
## 18 -74.29120 41.22047 2018-04-13T19:09:43.811Z NA 19.906
## 19 -74.29069 41.22057 2018-04-13T19:10:25.000Z 185 23.000
## 20 -74.29070 41.22057 2018-04-13T19:10:28.000Z 184 16.000
## 21 -74.29071 41.22058 2018-04-13T19:10:30.000Z 183 15.000
## 22 -74.29115 41.22049 2018-04-13T19:10:07.017Z NA 22.561
## 23 -74.29082 41.22050 2018-04-13T19:20:28.000Z 173 32.000
## 24 -74.29077 41.22047 2018-04-13T19:20:11.969Z NA 21.443
## 25 -74.29078 41.22048 2018-04-13T19:20:26.408Z NA 22.295
# Enter the row number of the point feature below
PointFeat<- 1
# Enter the row numbers of the line features below, separated with a colon
LineFeat<- 2:10
# As above, enter the row numbers of the polygon features below, separated with a colon
PolyFeat<- 11:22
# Now, use these row numbers to create individual 'spatial' objects for each of the features
# Point features
PointShape <- SpatialPointsDataFrame(wps[,1:2][PointFeat,], data = wps[PointFeat,])
# Line features
PointLine <- SpatialLinesDataFrame(sl = SpatialLines(list(Lines(list(Line(coords = wps[,1:2][LineFeat,])), ID='a'))), match.ID = F, data = as.data.frame("Line"))
# Polygon features
PointPoly <- SpatialPolygonsDataFrame(Sr = SpatialPolygons(list(Polygons(list(Polygon(wps[PolyFeat,][,1:2])),"Poly"))), data = data.frame("Poly"), match.ID = F)
Saving as ESRI shapefiles
writeOGR(obj = PointShape, dsn = paste(getwd()), layer = "PointFeature", driver = "ESRI Shapefile")
writeOGR(obj = PointLine, dsn = paste(getwd()), layer = "LineFeature", driver = "ESRI Shapefile")
writeOGR(obj = PointPoly, dsn = paste(getwd()), layer = "PolyFeature", driver = "ESRI Shapefile")
Test for accuracy
PointShape
ggmap(manmap) +
geom_point(data = data.frame(PointShape@coords),
color = "red",
size =1) # for better resolution, change size to smaller
# LineShape
ggmap(manmap) +
geom_polygon(data = PointLine,
aes(x = long, y = lat),
color = "red",
size =1) # for better resolution, change size to smaller
# PolyShape
ggmap(manmap) +
geom_polygon(data = PointPoly,
aes(x = long, y = lat),
color = "red",
size =1) # for better resolution, change size to smaller
## Regions defined for each Polygons
